Comparaison du MIC de plusieurs matrices de spectres

Comparaison du MIC de matrices de spectres Données human Serum prétraitées par SOAP par deux méthodes:

1. Full : Peps avec l’ensemble des étapes
2. Manuel : étapes appliquées manuellement avec Topspin par l’équipe ULG

1 Lecture des matrices des spectres et initiatilisation des variables utiles

#  Choix du répertoire où sont les matrices de données
pathdir="~/Dropbox/Partage_OMICS_COURSE/Data_stat2340/Human_Serum"
#  Definition des matrices à lire dans ce répertoire
data_used=c("HumSerum_full_PEPS.csv","HumSerum_Manual_500P.csv")
spectral_matrices_names=c("HumSerum_full_PEPS","HumSerum_Manual_500P")
nmatrices=length(data_used)
# Lecture des matrices de données.  
## On les stocke dans spectral_matrices_list
## On remplace les noms des colonnes par des noms sans le X
spectral_matrices_list=list() 
for(i in 1:nmatrices)
{
spectral_matrices_list[[spectral_matrices_names[i]]]=as.matrix(read.table(file.path(pathdir,data_used[i]),header=T,sep=";",row.names=1,dec="."))
dimnames(spectral_matrices_list[[i]])[[2]]=as.numeric(substring(dimnames(spectral_matrices_list[[i]])[[2]],first=2))
cat("\n Lecture de la matrice : ",spectral_matrices_names[i])
}

## 
##  Lecture de la matrice :  HumSerum_full_PEPS
##  Lecture de la matrice :  HumSerum_Manual_500P

# On prend la matrice des spectres communs dans les 2 bases (si ce n'est pas le cas au départ), on eleve les blancs des noms 
dimnames(spectral_matrices_list[[1]])[[1]]=gsub(" ","",dimnames(spectral_matrices_list[[1]])[[1]]) # on enleve des blancs des noms des spectres.  
namespectra=dimnames(spectral_matrices_list[[1]])[[1]]
countspectra=matrix(nrow=nmatrices+1,ncol=2)
dimnames(countspectra)=list(c(spectral_matrices_names,"Common Spectra"),c("Nb of spectra","Nb Buckets"))
countspectra[1,1]=length(namespectra)
countspectra[1,2]=dim(spectral_matrices_list[[1]])[2]
for(i in 2:nmatrices)
{
dimnames(spectral_matrices_list[[i]])[[1]]=gsub(" ","",dimnames(spectral_matrices_list[[i]])[[1]]) 
newnames=dimnames(spectral_matrices_list[[i]])[[1]]
compnames=is.element(namespectra,newnames) 
countspectra[i,1]=length(newnames)
countspectra[i,2]=dim(spectral_matrices_list[[i]])[2]
namespectra=namespectra[compnames]
}  
countspectra[nmatrices+1,1]=length(namespectra)
pander(countspectra)

	Nb of spectra	Nb Buckets
HumSerum_full_PEPS	32	500
HumSerum_Manual_500P	31	500
Common Spectra	31	NA

# On extrait les bons spectres et on les met dans le même ordre
for(i in 1:nmatrices)
{
spectral_matrices_list[[i]]=spectral_matrices_list[[i]][namespectra,]
} 
# Creation des variables qui définissent les groupes
# Pour les données human serum c'est le 5ième charactère du nom du spectre
groupes=as.numeric(substring(dimnames(spectral_matrices_list[[1]])[[1]], 5, 5))
names_spectra=substring(dimnames(spectral_matrices_list[[1]])[[1]],1,5)
repetition=as.numeric(substring(dimnames(spectral_matrices_list[[1]])[[1]], 2, 2))

# 3. on dessine les spectres
for(i in 1:nmatrices)
{
print(spectral_matrices_names[i])  
Draw(spectral_matrices_list[[i]],type.draw="signal",subtype="stacked",num.stacked=4)
}

## [1] "HumSerum_full_PEPS"

## [1] "HumSerum_Manual_500P"

2 Etude de la répétabilité spectrale

2.1 Between and Within group inertia

# Impression des résultats des calculs d'inertie
cat("Inerties non standardisées")

## Inerties non standardisées

pander(list(Inerties_non_standardisees=Inertiavalues,Inerties_en_PC=Inertiapc,Inerties_within_groups=InertiaPerGroup))

• Inerties_non_standardisees:

  	BI	WI	TI
  HumSerum_full_PEPS	6366	702.3	7069
  HumSerum_Manual_500P	0.02525	0.09521	0.1205

• Inerties_en_PC:

  	BI	WI	TI
  HumSerum_full_PEPS	90.06	9.936	100
  HumSerum_Manual_500P	20.96	79.04	100

• Inerties_within_groups:

  	HumSerum_full_PEPS	HumSerum_Manual_500P
  1	168.8	0.02751
  2	256.4	0.0251
  3	84.12	0.02492
  4	192.9	0.01768
  Total	702.3	0.09521

2.2 PCA

2.2.1 Eigenvalues

2.2.2 PCA Score Plots

draw.scores12 = vector("list")
draw.scores34 = vector("list")
for(i in 1:nmatrices)
{
draw.scores12[[spectral_matrices_names[i]]]=MBXUCL::DrawScores(PCA.res[[i]], drawNames=TRUE, type.obj = "PCA",createWindow=FALSE, main = paste0("PCA score plot for ",spectral_matrices_names[i]),pch = groupes, col = groupes,axes =c(1,2))
draw.scores34[[spectral_matrices_names[i]]]=MBXUCL::DrawScores(PCA.res[[i]], drawNames=TRUE, type.obj = "PCA",createWindow=FALSE, main = paste0("PCA score plot for ",spectral_matrices_names[i]),pch = groupes,col = groupes, axes =c(3,4))
}

draw.scores12

## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

draw.scores34

## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

2.2.3 PCA Loading Plots

draw.loadings12 = vector("list")
draw.loadings34 = vector("list")
for(i in 1:nmatrices)
{
draw.loadings12[[spectral_matrices_names[i]]] = MBXUCL::DrawLoadings(PCA.res[[i]],  type.obj = "PCA",createWindow=FALSE, main = paste0("PCA loadings plot for", spectral_matrices_names[i]),axes = c(1:2), loadingstype="l", num.stacked = 2)[[1]]
draw.loadings34[[spectral_matrices_names[i]]] = MBXUCL::DrawLoadings(PCA.res[[i]],  type.obj = "PCA",createWindow=FALSE, main = paste0("PCA loadings plot for", spectral_matrices_names[i]),axes = c(3:4), loadingstype="l", num.stacked = 2)[[1]]
}
draw.loadings12

## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

draw.loadings34

## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

2.3 Unsupervised clustering (MIC)

nClust = length(unique(groupes)) # nombre de clusters à rechercher
clustres=matrix(0,nrow=nmatrices,ncol=8)
dimnames(clustres)=list(spectral_matrices_names,c("DunnW","DunnKM","DBW","DBKM","RandW","RandKM","AdjRandW","AdjRandKM"))
for(i in 1:nmatrices)
{
print(spectral_matrices_names[i]) 
ClustMIC.res = MBXUCL::ClustMIC(Intensities = spectral_matrices_list[[i]], nClust = nClust, Trcl = groupes, Dendr = TRUE)
clustres[i,]=as.numeric(ClustMIC.res[1,1:8])
}

## [1] "HumSerum_full_PEPS"

## [1] "HumSerum_Manual_500P"

pander(clustres)

Table continues below

	DunnW	DunnKM	DBW	DBKM	RandW	RandKM
HumSerum_full_PEPS	0.8246	0.8246	0.6776	0.6776	1	1
HumSerum_Manual_500P	0.9332	0.9332	0.7824	0.7824	0.7054	0.7054

	AdjRandW	AdjRandKM
HumSerum_full_PEPS	1	1
HumSerum_Manual_500P	0.2324	0.2324

2.4 PLS-DA

nLVPLSDA = 4  # nombre du composantes du PLSDA
nrep=nlevels(as.factor(groupes)) # Nombre de réponses
PLSDA.res=list()
perf.plsda.RMSEP=matrix(nrow=nmatrices,ncol=nrep)
dimnames(perf.plsda.RMSEP)=list(spectral_matrices_names,paste0("Y",1:nLVPLSDA))
perf.plsda.R2=perf.plsda.RMSEP
perf.plsda.Q2=perf.plsda.RMSEP
for(i in 1:nmatrices)
{
cat("PLS-DA pour ",spectral_matrices_names[i])  
PLSDA.res [[i]]= PLSDA(x = spectral_matrices_list[[i]], y = groupes, nLV = nLVPLSDA, drawRMSEP = TRUE)
perf.plsda.RMSEP[i,] = PLSDA.res[[i]]$RMSEP
perf.plsda.R2[i,]=PLSDA.res[[i]]$R2
perf.plsda.Q2[i,]=PLSDA.res[[i]]$Q2
}

## PLS-DA pour  HumSerum_full_PEPS

## PLS-DA pour  HumSerum_Manual_500P

pander(list(RMSEP=perf.plsda.RMSEP,R2=perf.plsda.R2,Q2=perf.plsda.Q2))

• RMSEP:

  	Y1	Y2	Y3	Y4
  HumSerum_full_PEPS	0.1145	0.0792	0.1061	0.1117
  HumSerum_Manual_500P	0.3953	0.2478	0.426	0.2985

• R2:

  	Y1	Y2	Y3	Y4
  HumSerum_full_PEPS	0.9728	0.9873	0.971	0.9662
  HumSerum_Manual_500P	0.619	0.8414	0.6475	0.6968

• Q2:

  	Y1	Y2	Y3	Y4
  HumSerum_full_PEPS	0.9316	0.9672	0.9356	0.9349
  HumSerum_Manual_500P	0.1839	0.6794	-0.0381	0.5347

2.4.1 Score plots de la PLS_DA

PLSDA.scores= vector("list")
for(i in 1:nmatrices)
{ 
PLSDA.scores[[spectral_matrices_names[i]]]=DrawScores(PLSDA.res[[i]], drawNames=TRUE, type.obj = c("PLSDA"),
        createWindow=FALSE, main = paste0("PLSDA score plot for ", spectral_matrices_names[i]),pch = groupes,col = groupes, axes =c(1,2))
}
PLSDA.scores

## $HumSerum_full_PEPS

## 
## $HumSerum_Manual_500P

2.4.2 Loading plots de la PLS_DA

PLSDA.loadings= vector("list")
for(i in 1:nmatrices)
{ 
PLSDA.loadings[[spectral_matrices_names[i]]]=DrawLoadings(PLSDA.res[[i]],  type.obj = "PLSDA",
        createWindow=FALSE, main = paste0("PLSDA loadings plot for", spectral_matrices_names[i]),
        axes = c(1:2),  loadingstype="l", num.stacked = 2)
}
PLSDA.loadings

## $HumSerum_full_PEPS
## $HumSerum_full_PEPS[[1]]

## 
## 
## $HumSerum_Manual_500P
## $HumSerum_Manual_500P[[1]]